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Abstract 

We describe a constructive procedure to separate overlapping infrared divergences in multi- 
loop integrals. Working with a parametric representation in _D = 4 — 2e dimensions, adequate 
subtractions lead to a Laurent series in e, where the coefficients of the pole- and finite terms 
are sums of regular parameter integrals which can be evaluated numerically. We fully autom- 
atized this algorithm by implementing it into algebraic manipulation programs and applied it 
to calculate numerically some nontrivial 2-loop 4-point and 3-loop 3-point Feynman diagrams. 
Finally, we discuss the applicability of our method to phenomenologically relevant multi-loop 
calculations such as the NNLO QCD corrections for e^e^ —> 3 jets. 
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1 Introduction 



The increasing experimental precision at present and future colliders requires fast progress in the 
calculation of higher order corrections from the theoretical side. A crucial role is thereby played 
by the calculation of multi-loop integrals, which becomes an increasingly challenging task as the 
number of loops and the number of kinematic invariants gets larger. 

In the case of two-point functions, up to four loop orders could be evaluated ^ by exploiting 
the integration-by-parts method Q and recurrence relations. This powerful technique also has been 
generalized to calculate three-loop propagator- type diagrams in heavy quark effective theory Q. 
Similarly, the presence of only one kinematic invariant allowed for the calculation of two-loop three- 
point functions [|[ ||, 0, |[ which could be applied to compute the partonic cross sections 
for DIS and the Drell-Yan process to NNLO |ll|. Physical processes with a more difficult 
kinematic structure could not be computed to two-loop order yet, amplitudes exist only for very 
special cases [Q. A major reason is that two-loop integrals with four external legs, involving more 
than one kinematic variable, constitute a more complex problem. Only very recently, the analytic 
calculation of the massless planar and non-planar double box could be achieved by using 
Mellin-Barnes integration techniques. Furthermore, results for single-box integrals with self-energy- 
and vertex insertions have been given in and techniques based on differential equations to reduce 
two- loop four point functions to a set of master integrals have been derived in |l^ . The tensor 
reduction of massless double box integrals has been completed in where the second master 
integral needed in the reduction of crossed two-loop boxes is given. However, it is not clear at the 
moment whether the techniques developed so far are sufficient to solve more complicated higher 
order problems as well. 

Hence, to obtain results for integrals which are beyond the scope of analytic integration methods, 
but also as a simple check of involved analytic calculations, it is desirable to have a method at hand 
which allows the calculation of multi-loop integrals numerically. For two-loop box integrals with 
internal masses, thus with no infrared divergence, semi-numerical approaches have been designed 
in The massless case requires a different approach because of the presence of infrared diver- 

gences, which first have to be extracted from the integral in order to make it amenable to numerical 
evaluation. 

In perturbative QCD, the calculation of infrared safe quantities has to be organized such that the 
infrared poles stemming from virtual and real higher order corrections cancel. At next-to-leading 
order usually semi-analytical methods are used to achieve this cancelation. As an alternative, a 
completely numerical method has been developed in pO| . The algorithm allows for the summation 
of the contributions from different cuts of a graph before integration, such that the cancelation of 
soft and coUinear divergences is built in as imposed by unitarity. Nevertheless, since the method 
requires the deformation of the multi-dimensional integration contours to avoid fake singularities, 
its general applicability to NNLO calculations is highly nontrivial. 

Hence a general local subtraction procedure to separate infrared divergences from individual 
graphs of arbitrary loop order would be very useful. The IR singularity structure of higher-loop 



Feynman diagrams was investigated in pi[ in four dimensions, and later in |22 in the context 



of dimensional regularization and factorization in QCD. Working in dimensional regularization, 
subtraction procedures are well known for UV poles, and also for IR divergences present in Euclidean 
space (see e.g. |2^]). On the other hand, no general subtraction scheme for soft and coUinear IR 
singularities arising in Minkowski space is known for individual graphs. The method we present in 
this paper has been designed to isolate poles in the dimensional regulator e for an arbitrary Feynman 
graph. Although the method also works for one-loop integrals, its virtues show up rather in two- or 
higher loop integrals with iV > 3 external legs, at least one of them being massless. It allows one to 
disentangle overlapping soft and coUinear divergent regions in Feynman parameter space by dividing 
the latter into sectors where parameters can get singular only in an independent manner. Then, 
by adding and subtracting adequate counterterms, one can isolate the singular parts and perform 



1 



the integrations over the corresponding parameters analytically. The remaining regular integrals 
are in general too complex for analytical integration, but they can be integrated numerically. This 
procedure is quite general and can in principle be applied to graphs with an arbitrary number of 
loops and legs, the limitations being only disk space and computing time. 

The algorithm can be divided into four blocks. In the first one, the J-distribution constraint 
on the Feynman parameters is eliminated in a particular way and in the second one, the singular 
contributions are isolated in parameter space. The third block consists of a subtraction procedure 
for the 1/e™ poles, producing a set Cm of finite functions of the Feynman parameters as coefBcients 
of each order-m pole. In the fourth block, the integrations over the Feynman parameters in these 
functions arc performed. The coefficient functions of the leading and sublcading pole of a graph 
are in general simple enough to be integrated analytically with an algebraic manipulation program, 
whereas the remaining functions have to be integrated numerically. 

The paper is organized as follows. In section two. wc outline the algorithm. In section three, 
some examples are given in order to show the applicability of the method to quite different types 
of integrals. First we treat some two-loop four-point functions with one external leg off-shell as 
an example for a three-scale problem which could not be solved analytically yet. Then three- loop 
three-point graphs with two external legs on-shell and eight and nine propagators, respectively, are 
calculated. Section four contains a discussion of available analytical/numerical results for elements 
entering the c;ak;idation of various processes of phenomenological relevance, pointing out where our 
method improves the situation. 



2 The algorithm 



In this section we describe the algorithm to treat a general £)-dimensional scalar L loop Feynman 
diagram with N propagators. If E is the number of external legs with momenta {pi, . . .ps}, L the 
number of loop momenta {fci, . . . , Aii,}, {mi, . . .ttin} the (not necessarily nonzero) masses of the 
propagators Pje{i,...,N}j and dlCjne{i,...,L} = d^km/[iT^''^^^^, a general scalar Feynman diagram G 
can be written as 



N 

G = / d/Ci . . . d/Ci [] Pj{{k}, {p}, m^) 



(1) 



In the present paper we will not deal with powers of propagators different from one. Nevertheless, 
higher (and even non-integer) powers of propagators can be treated with basically the same algo- 
rithm. Introducing Feynman parameters, the integral can be expressed in terms of a symmetric 
{L X L)-matrix M, an L-vector Q (with 4-vectors in each component) and a scalar function J. The 
contraction of Lorentz indices is indicated by a dot. 

n -AT 

L L 

Y,kj-kiMji-2j2kj-Qj+J (2) 



G = r(iV) j d^x S{l-^xi) j dJCi-.-dJCL 



After having shifted the loop momenta to get rid of the linear term, one obtains after momentum 
integration the following parameter representation of the graph G. 



G = {-lfT{N-LD/2) d^xS{l-J2xi) 



N 



UN-{L+l)D/2 
jrN-LD/2 



(3) 



where 



J^{x) = det(M) 



3,1=1 



(4) 
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U{x) = det(M) 



(5) 



As is well known, the parametric representation is determined by two functions which we call U and 
!F. They are defined by the topology of the corresponding Feynman diagram |2j, ^ . contains 
the Mandelstam variables related to the different cuts of the graph. 

A necessary condition for the presence of infrared divergences is that the Landau equations ^ 
are fulfilled. A representation of the Landau equations following from Eq. (0) is given by 



T = (6) 



The function lA cannot lead to infrared divergences of the graph, since giving a mass to all external 
legs would not change U. Apart from the fact that the graph may have an overall UV divergence 
contained in the overall F-function (see Eq. @), UV subdivergences may also be present. A necessary 
condition for these is the vanishing of U. In this way the UV poles can be identified and treated 
according to standard renormalization procedures. Infrared poles in 1/e come from the parameter 
region where some Feynman parameters are small such that J- vanishes. The Feynman parameter 
integrations which lead to poles can be related to a kinematical configuration in momentum space 
by using Eqs. (^. The IR poles come from soft and/or coUinear momentum configurations where 
propagators do not describe virtual particles anymore but rather the propagation of an on-shell 
particle together with eventual splittings^. The soft /collinear poles of multi-loop graphs are the 
result of such kinematical situations in loop momentum space. These singular regions are generally 
not separated from each other in momentum (or equivalently in Feynman parameter space), they 
are overlapping. 

Our aim is to disentangle these regions of overlapping IR divergences. To this end we use a 
method called sector decomposition^ Iterated application will lead to a set of integrals in which 
the infrared singular behaviour is not contained in complicated functions anymore, but in simple 
products of Feynman parameters raised to some power, times remnants of the functions J-', lA, whose 
structure is such that they neither lead to a pole anymore, nor change the exponent of the respective 
poles. In more detail, the procedure consists of four basic building blocks: 

Part I Generation of primary sectors 

In the first part of the algorithm we split the integration domain into N parts and eliminate the 5- 
distribution in such a way that the remaining integrations are from to 1. To this end we decompose 
the integration range as follows 



N N ,ryn N 



/ d^x= d^x\[e{xj>{))^^ d^xWe{xi>Xj> 

Jo Jo ,_^ Jo . , 



0) (7) 



j=l 1=1 " " J=l 



where the ^-function is defined as 

01 ^ , _ r 1 if a; > y is true 
~ \ otherwise 

The integral is now split into N domains corresponding to N integrals Gi from which we extract a 
common factor: G = (— 1)^F(A^ — LD/2) J2iLi Gi- In the integrals Gi we integrate out xi by using 

^This can be visualized by defining the reduced graph of a diagram, which is the diagrammatic representation 
of solutions of the Landau equations. 

^This method was of some importance in the history of UV regularization, i.e. to establish the BPHZ method. It 
was used by Hepp [ p9| to deal with overlapping UV divergences. 
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the (5-distribution after having done the substitution 




, j <l 

, J=l (8) 
, j>l 

Because of homogeneity, Xi factorizes completely in the functions V({x) — > Ui{t)x\' and J^{x) 
Ti{t)xf~^^ and thus, using J dxi/xi5{\ — xi{\ + X]fe=L^*fe)) = 1' obtains 

1 j.N-iL + l)D/2 

Gi = j d'^-'t^^W^LDj^ , 1^1,. ..N (9) 
' 

Note that the singular behaviour leading to e-poles still comes from the regions of small i's. This 
feature would be lost if one integrated out the i5-distribution in a naive way, since this would produce 
poles at upper limits of the parameter integral as well. The generated sectors will be called primary 
sectors in the following. The functions Ui and J-i are polynomials in the parameters tj . 

Part II Iterated sector decomposition 

The second part of the algorithm consists of the iterated application of sector decomposition and 
a remapping of parameter space to the unit cube in order to disentangle the overlapping singular 
regions of the integrands. Starting with Eq. (^) one repeats the following steps until complete 
separation of overlapping regions is achieved. 

II. 1: Determine a minimal set of parameters, say S = {tai, • ■ • ,tar}: such that lAi, respectively J-i, 
vanish if the parameters of S are set to zero. S is generally not unique. Additional selection 
criteria can be introduced to choose an S which does not lead to a large number of subsequent 
sector decompositions. 

II. 2: Decompose the corresponding r-cube into r subsectors. 

r r r 

n > > 0) = 5] n ^(^"'^ ^ ^ 0) (10) 

j = l fc=l i=i 

II. 3: Remap the variables to the unit cube in each new subsector by substituting 

^ J, J ^Qfc^Qj J j 7^ k (11) 

This gives a Jacobian factor of t^^. By construction ta^ factorizes at least from one of the 
functions Ui, Ti. The resulting subsector integrals have the general form 

\ (n-I \ N-{L + l)D/2 

Vj=i / -^ik 

For each subsector the above steps have to be repeated as long as a set S can be found such that 
one of the functions Ui,,,, Ti,,, vanishes if the elements of S are set to zero. In each subsector new 
subsectors are created, resulting in a tree-like structure after a certain number of iterations. The 
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book-keeping can be done with respective multi- indices. The iteration stops if the functions Uikik2...i 
^ikik2... contain a constant term, i.e. if they are of the foUowing schematic form 

Uik,k2... = l + u{t) (13) 

where u(t) and ff3{t) are polynomials in the variables tj (without a constant term), and S/j are 
kinematic invariants, defined through (^). Thus, after a certain number of iterations, each integral 
Gi is split into a certain number, say R, of subsector integrals. For simplicity we replace the multi- 
index fcifc2 . . . stemming from the subsector decomposition by a single index which just counts the 
number of generated subsectors. Now, the produced subsector integrals are exactly of the same 
form as in Eq. (p^, with the difference that the index k now runs from 1 to R, the total number of 
produced subsectors. 

Evidently the singular behaviour of the integrand now can be trivially read off the exponents 
Aj, Bj for a given subsector integral (Aj, Bj are integers). The singular behaviour is manifestly 
non-overlapping now and thus it is straightforward to define subtractions. Before doing so a few 
comments are in order. 

The described method cannot always lead to an optimal sector decomposition in the sense that 
the integral one starts with is split into the minimal number of subsector integrals, since obviously 
even finite integrals would be decomposed if a set S exists such that U or J- vanish. The virtue of the 
algorithm lies in its easy programmability, as we introduced standardized representations suitable 
for iteration. 

For massless 2-loop 4-point functions with 7 propagators and all external legs on-shell the number 
of generated subsectors is a few hundred. For 3-loop 3-point functions with two legs on-shell and 
9 propagators it is a few thousand. Hence the bookkeeping of such numbers is only possible on a 
computer. 



Part III Extraction of the poles 

The third part of the algorithm consists in a pole subtraction procedure for the subsector integrals 
Gik- As the infrared sensitive variables now factorize in the subsector integrals, one can separate 
the part of the integrand which leads to e-poles. 

Explicitly, the following procedure has to be worked through for each variable tj^i ^v-i and 

each subsector integrand: 

• The integrand of Eq. (|l^), characterized by the respective exponents Aj—Bje {j = 1, . . . , iV— 1) 
of tj and the functions of the form (|l^), can for each tj be written as 

I,= I dt,tf^-'''''>I{t,,e) (14) 



If Aj > 0, the integration does not lead to an e-pole. In this case no subtraction is needed 
and one can go to the next variable ij+i- If Aj < 0, one expands I(ij, e) into a Taylor series 

around tj = 0. Using the definition ^^^^(O, e) dPI{tj,e)/dtP 



, one obtains 

ti=0 



I{t„e)^ J2 if (0,6)^ + i?(t„6) (15) 



pi 
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• Now the pole part can be extracted easily, and one obtains 



p=0 

By construction the integral containing the remainder term R(tj,e) does not get poles in e 
from the tj -integration anymore. For example, in the generic case of a logarithmic divergence, 
Aj — —1, p — and R{tj, e) — 2{tj, e) — X,(0, e). Since, as long as j < — 1, the expression 
( p^ ) still contains an overall factor f^^^ ■^j+i'^^ jg ^j-^g g^j-Qg form as (|lj) for j — > j + 1 
and the same steps as above can be applied to it. 

After iV — 1 steps all singular integrations are done analytically and all poles are extracted. The 
resulting expression can be expanded in e now. This defines a Laurent series in e with coefficients 
Cik,m for each subsector integral Gik- Since each loop can contribute at most one soft and coUinear 
term, the highest possible infrared pole of an L— loop graph G is l/e^^. 

Gi>^-i:^+o{e) (17) 

m=0 

Following the steps outlined above one has generated a regular integral representation of the coeffi- 
cients Cik^m, consisting of {N — 1 — m)-dimensional finite integrals over parameters t. 

Symmetries of the graph typically lead to equalities between primary sectors. It is thus useful 
to calculate the primary sectors G/ 

R 

Gi = Y,Gik (18) 

k=a 

separately before summing over the I subsectors. In this way the symmetry relations provide a 
nontrivial check of the calculation. 



Part IV Calculation of the pole coefficients 

Part four of the algorithm consists in the computation of the finite subsector integrals. The integrals 
contributing to the leading pole give ratios of polynomials in the Mandelstam variables. For the 
coefficient of the subleading pole one generally gets logarithmic terms. In principle, one can attempt 
to perform the (N — 1 — r7i)-dimensional integrations of all functions contributing to the coefficient 
of the 1/e™ pole analytically. However, the sector decomposition produces many surface terms 
which will cancel only after summing up all subsector integrals, such that the analytical integrations 
become more and more involved for smaller values of m, especially if the graphs contain more than 
one scale, as for example in the case of 2-loop 4-point functions. For these functions only the 
coefficients of the leading and subleading poles could be obtained analytically by automatizing the 
integrations using Mathematica jss). Pushing the analytical integrations further is possible only if 
some of the more complicated functions are manipulated "by hand" before feeding them into the 
subroutine, but this is tedious in view of the large number of functions to integrate. More powerful 
analytical integration routines, specialized to manipulations of poly logarithms and Nielsen functions 
would be needed to allow for a complete analytical treatment. 

On the other hand, the parametric integral representations are all very well suited for numerical 
integration, as long as the parametric function has a definite sign. Then, it contains at most 
integrable, logarithmic divergences at the border of the integration domain. These generally present 
no problems for the numerical integrators which are on the market. If is not of a definite sign, 
which means that one has to integrate over thresholds, the integrands contain poles inside the multi- 
dimensional integration domain. The presence of these poles typically considerably slows down, if 
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not hinders at all, the numerical evaluation of the integral. More advanced integration algorithms 
like e.g. the one proposed in ||3^] may solve this problem. We restrict ourselves to the case of definite 
sign here. 

In most algebraic programs it is directly possible to create FORTRAN functions from a given 
expression. We fully automatized the translation of the expressions for the subsector integrals into 
the FORTRAN codes. For every primary sector we calculated the integral of the sum of the subsector 
integrands ( p^ with the Monte Carlo program BASES |3^. The integrations of all examples we 
calculated were totally stable. The limitation for the numerical integration comes thus only from 
CPU time, which increases with the number of loops and legs or, correspondingly, with the number of 
subsector functions. We will give more detailed information on program parameters in the examples 
below. 



3 Applications 

In this section we apply our method to calculate several nontrivial Feynman diagrams. To this 
end we implemented the algorithm outlined in section 2 into algebraic manipulation programs. To 
crosscheck the output we created two independent codes, one written in Maple the other one in 
Mathematica |Q . First we show a comparison of the results obtained by our method with the results 
from the analytical calculation of the planar jisj and non-planar iQ massless double box. Then 
we give results for some 2~loop 4-point functions for which analytical results are not yet available 
in the literature, relevant for the calculation of certain higher order QCD corrections (see section 
These are the massless planar and non-planar double box with one external leg off-shell, given 
in Figs. 1-3. As discussed above, our method so far only allows us to calculate numerical values of 
these graphs if all the Mandelstam variables have the same sign. In any case our result may serve 
as a nontrivial check of a future analytical computation of these graphs. Further we will give results 
for 3-loop 3-point graphs with two on-shell legs. 



The numerical values given in subsection 3.1 contain a relative error of one percent. Note 



that for comparison purposes our conventions for the prefactors (F-functions) should be used since 
multiplying our results with conversion factors that are a power series in e may lead to a bad error 
propagation^. In subsections 3.2 and 3.3 we use a different prefactor that in subsection 3T because it 



leads to more compact expressions for the analytical result. This conversion may lead to errors in the 
numerical result which are slightly larger than one percent due to the error propagation mentioned 
above. 

The numerical calculations were done on a DEC-ALPHA workstation running with an EV6 
processor. The CPU time needed ranges from about 5 hours for the planar two-loop graph up 
to about 3 days for the 3-loop graph with 9 propagators. The computer time needed is always 
dominated by the finite 0{e'~') terms. 



3.1 Comparison with analytical results 

In this subsection we compare our numerical approach with the analytical results for massless two- 
loop box diagrams calculated only recently ||l^, These results have also been crosschecked 
analytically meanwhile, and the perfect agreement with our results within the error of numerical 
integration confirms the reliability of our method. 



The massless planar double box Bj{s,t) 

The graph (see Fig. |l|with leg four on-shell) depends on the two kinematical invariants s = (pi+P2)^ 
and t = {p2 + P^Y ■ We show a comparison to the analytical result at two numerical points, the 

^ The reason is that multiplication with such a conversion factor mixes the different pole coefficients in our Laurent 
series, which can lead to a larger relative error in the converted result. 
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symmetric point s = t = —1 and the asymmetric point {s,t) = (—1,-1/2). As stated above, the 
physical situation st < is not suited for numerical evaluation in a straightforward manner because 
thresholds are present. In the following table we show the comparison between the analytical result 
of ||l3| and our numerical result for the relevant Laurent coefficients c^(s,t) of the massless planar 
double box Bf{s,t), defined through Eq. (^ and the subtraction algorithm: 

Bf{s,t) = (-iyr{3 + 2e)J2^^^ + Oie) (19) 

m=0 

One can check that the agreement is always better than the demanded one percent. 









-1,-1) 








C4 


cr 




Ci 




analytical 


-2. 


6. 


4.9167 


-11.495 


-13.801 


numerical 


-2.0000 


6.0000 


4.9188 


-11.492 


-13.811 








-1,-2) 








C4 




^2 


Ci 




analytical 


-1. 


3.8664 


-0.38116 


-9.2384 


-2.9973 


numerical 


-1.0000 


3.8664 


-0.38059 


-9.2377 


-2.9990 



Table 1: The massless planar double box. Comparison between analytical and numerical result for 
the points {s,t) = (—1,-1) and (—1,-2). 



The massless non-planar double box B^^{s,t,u) 

Without using momentum conservation, the graph (see Fig. ^ with leg one on-shell) depends on the 
three kinematical invariants s — {pi t — {p2 +^3)^ and u = {pi +pzY ■ In the physical region 

one invariant is positive, leading always to an imaginary part. To avoid the corresponding threshold 
in the numerical integration we treated s,t,u as independent parameters. We compared to the 



analytical result of |14| by calculating numerically the coefficients c'^^{s,t,u) of the Laurent series 
of Bj^{s, t, u) (defined with the same conventions as in (p^)) at the symmetric point s = t = u = — 1 
and the asymmetric point (s, t, u) = (—1, —2, —3), see Table 2. 



3.2 Two— loop massless 4— point functions, one leg off-shell 

Now we turn to the calculation of graphs where fully analytical results do not exist yet. The 
Mandelstam variables s, t, u are defined as above. If an external leg pk is off-shell or on-shell but 
massive, we write = m^. For the 4-point functions under consideration, exactly one leg is off- 
shell. Then momentum conservation implies s + i + u — mf,, but it has to be emphasized that this 
constraint has not been used to obtain the analytical results, in order to be able to compare them 
to numerical results for unphysical kinematics as well. 

For sums of Feynman parameters we use short-hand notations like e.g. Xi -f Xj + Xk+xi — Xijki- 
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{s,t,u) = (- 


1,-1,-1 


) 






C4 




^2 




„NP 



analytical 


-1.75 


3. 


22.828 


-113.63 


395.26 


numerical 


-1.7500 


2.9960 


22.818 


-113.75 


393.08 




{s,t,u) = (- 


1,-2,-3) 






C4 








Co 


analytical 


-0.4167 


0.9310 


5.8586 


-42.760 


162.81 


numerical 


-0.4167 


0.9295 


5.8748 


-42.614 


164.16 



Table 2: The massless non-planar double box. Comparison between analytical and numerical result 
for the points (s, t,u) = (—1, —1, —1) and (—1, —2, —3). 
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Figure 1: The planar double box with leg 4 off-shell. 

The graph t, u, m|) 

With the labeling as in Fig. |l| one finds for the functions W, J-: 

W = 2^1232:567 + XiXi2S,^(,7 

T = i-s){x2X3X45er + X5XeXi234: + X2X4,X6 +X3X4X5) 

+ {-t)xiX4^X7 + {-■m'l)x-j{x2X^ + X;^Xi23i) (20) 

The sector decomposition produces about 200 subsector integrals. For the leading and subleading 
pole we get the following analytical result: 

Bl.mass = r^(l + e)(-m^)-2'-i^(^l-|[log(sM)+log(Vm^)])+0(l) (21) 

Numerically we find for the points (-1/3, -1/3, -1/3, -1) and (-1/2, -1/3, -1/6, -1): 

,„.„,(-l/3. -1/3. -1/3. -1) . r»(l + (-™ - - - 2^ - 1«,) ,22) 

_H:000 _«^_58|8_2a85^ ^ 

e'* e-' e / 

The graph B^f^^,,^^{s,t,u,ml) 

With the labeling as in Fig. |^, the functions U, T are given by: 

U = 2:123X4567 + 2:45X57 



1 Amass \ 



-1) = r2(i 
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Figure 2: The non-planar double box with leg 1 off-shell. 



+ {—t)xiX4X7 + {~u)xiX5X(i + {-ml)xi{xQXT + X3a;4567) 



(24) 



We note that this graph contains hnear IR divergences, as it is the case if all legs are massless [Q. In 
a gauge theory, numerator functions would be present to prevent singularities of such strength. Since 
the subtractions are more complicated if linear divergences are present, the subsector integrands are 
more complicated as well. This leads to larger FORTRAN functions and thus slows down the 
numerical computation in comparison with the similar graph B^f^^^^ ^ below, which has poles 
stemming only from logarithmic divergences. It also requires a refined analytic integration routine. 
The sector decomposition produces about 250 subsector integrals. For the leading and subleading 
pole we obtain the following analytical result: 



nNl 

^7,1 



NP 

mass, a 



1 r 
73 



As'^tu ^ e"* 
'i[s + t + u + m\) 



s + i + 7i — m, H 7j + - [t + U) 

mi 2 



tu 



+ log(s/mf ) [-2(s + t + u-m{)- 2^ - {t + u)] 



t u 



+ log(t/mJ) [-2(s + t + u-ml) + 2— + it - u] 



1 

tVL 1 ^ 

+ log(M/rn?) [-2(s + t + u- m\) + 2^ + 3it - t] \ 

mi i ) 

Numerically we find for the points (-1/3, -1/3, -1/3, -1) and (-1/2, -1/3, -1/6, -1): 



T.lmass.a V 



r^(i 



5.997 101.7 592.7 3340. 



18522, 



(25) 



(26) 



rjJVi 
^7,1 



NP 

mass 



5.504 87.98 296.6 1753. 



11741, 



(27) 



The graph B^f^^^^,,{s,t,u,ml) 

With the labeling as in Fig. ^ one obtains for the functions U, T: 

U = a;i23a;4567 + XAf^x^^ 

T = i-s){x2X3X45e7 + X2Xexr -\- X3X4X5) 

+ { — t)xiX4X7 + { — u)xiX5Xq 

+ i-'ml){x,iX7Xi2346 + X2X4XJ + X^X^Xq) 



(28) 
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2 



Figure 3: The non-planar double box with leg 3 off-shell. 



The sector decomposition produces about 230 subsector integrals. For the leading and sublcading 
pole we obtain the following analytical result: 

Zstu ^ Ze^ 
+ ^\-{t + u) log{s/ml) 

+ (4m3 -t + u) \og{t/ml) 

+ {4ml + t-u)\og{u/ml)'\} + Oi^) (29) 

Numerically we find for the points (-1/3, -1/3, -1/3, -1) and (-1/2, -1/3, -1/6, -1): 
j^p 111 ,.-^2,^ ^( 31.50 108.8 2.560 779.^ 



57:r™a..,6(- 3, -3, -3, -1) = r^(l + e) ( ^ — + — + 2395. ) (30) 

r^NP /I 1 1 ,N J 40.50 203.9 357.3 409.9 

B^fmass.bi-^^ -3, - g, -1) = r'(l + e) -4 — + + 4608. j (31) 

We recall that all the numbers given were calculated for unphysical kinematics (all particles ingoing) 
in order to have positive definite denominators. The numerical results for the 1-mass two- loop graphs 
as given above do not correspond to a physical situation. Nevertheless, they can easily be related to 
a physical process by the following reasoning. In the case of a 1 — > 3 process, e.g. a virtual particle 
with a squared momentum of decaying into 3 massless particles, s,t,u and m? are positive. 
By factoring out {—m?) in the respective functions one gets again positive definite integrands. 
Especially, for a one-mass 2-loop graph with N propagators, one finds 

G2-'°°f (s, t, u, m^) = (-m^ - i5y^+^-^'G'^-'^°°P{s/m^ ,t/m^ ^u/m^ , 1) 

- 1 - - (-1)" ^'nG^~-^.^s/^\ "/™^ d 02) 

We introduced an infinitesimal imaginary part i5 where necessary. This shows that our method 
allows to calculate 2-loop scalar integrals which appear for example in the calculation of the NNLO 
QCD corrections of e+e^ 3 jets, as explained in more detail in section ^ 

3.3 Three— loop massless 3— point functions, one leg off-shell 

Now we want to present results for 3-loop 3-point integrals with two legs on-shell. These contain 
only one scale which can be factored out of the integrand, such that one has to calculate a pure 
number which can be done numerically once and forever. We restrict ourselves to only two examples, 
the graphs shown in Figs. ^ and ||. 
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Figure 4: A 3-loop Mercedes-Star (MS) topology. 



The graph M^gl^) 

With the labeling as in Fig. ^, the functions U, T read: 

U = X23478a;5a;i6 + a;234(a;6a;8 + a^ia^es + a;7a;568) + a;ia;7X68 + a;8(xia;6 + 2;5a;7) 
T = (-s)[xiX5(a;6a;23478 +a;7a;8) +a;iX4(a;6a:8 +X7a;568) 

+X2X5(a;8a:i67 + xqx-;) + a;2a;4(xi7a::58 + a;6a;i578)] (33) 

The sector decomposition produces about 1000 subsector integrals. We note that because of the 
symmetries of the graph only 5 of the 8 primary sectors as defined in Eq. have to be calculated, 
i.e. M5'8°'^i = MSl'"'^ MSl^"^ = MSl'"'^ MS^^"^ = MS^'"'^ On the other hand, the recalculation 
of these identical sectors are a good check of the algebraic/numerical computer routines. For the 
leading, subleading and subsubleading pole we find the following analytical result: 

Numerically we find: 



r3(l + e) /0.02778 0.0000 0.2288 0.6692 0.6152 2.005 ^„ , , 
MSs{s) = j-^^-^[^^^ + ^^ + — -^ + ^ + 17.85 ) (35) 



The graph BBTg^s) 




5 


2 








Pi 


4 


1 








P2 


6 


3 





Figure 5: The 3~loop box-box-triangle (BBT) graph 
With the labeling as in Fig. ^, the functions W, T are given by: 

U = a;i234a;7X89 + 0:789 (Xi23a;456 + 2:4X56) 

T = (-s) [2:20:3 (a;4567a;89 + 2:72:456) + (2:32:5 + 2:2X6)^42:789 + (x3a:8 + 2:22:9)^42:7 

+X5X6a:i2342:789 + (2:6X8 + X5a:9)a;7Xi234 + X8X9(a:i234X567 + X123X4)] 



(36) 
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The sector decomposition produces about 2400 subsector integrals. For the leading and sublcading 
pole we obtain the following analytical result: 

r^(i + .) ( 1 



BBTg{s) 



(-S - i6) 



3+3c 



36e6 



Numerically we find: 

r3(l + e) 



4 Discussion 



0.02778 0.0000 0.6852 2.072 6.613 25.07 



40.42 (37) 



We want to discuss now the phenomenological applicability of our method to standard QCD pro- 
cesses. We will focus on reactions with massless particles in the loop. The first step of a multi-loop 
calculation is to express the corresponding amplitudes, preferably decomposed in colour and helicity 
space, in terms of basic tensor and scalar integrals. In the analytical approach one tries to express 
all tensor integrals by a certain set of scalar integrals, so-called master integrals. In a numerical 
approach it may be more convenient to calculate certain tensor integrals directly. In a Feynman 
parameter language this amounts to the calculation of integrals of the type (^, where one has 
in addition a polynomial A of Feynman parameters in the numerator. The generalization of our 
method to include these nontrivial numerators is straightforward. One way to proceed is to carry 
along the corresponding polynomial A through all the steps in addition to J- and lA, iterating the 
decomposition until Aik contains also a constant term. This leads to an expression of type (12) with 
the integrand multiplied by Aik- Part III and IV then work in exactly the same way as above. Since 
the presence of numerators in general improves the IR behaviour, the number of necessary subsector 
decompositions will be reduced, which typically shortens the computation time of the diagram under 
consideration as compared to the scalar case. 



process 


kinematics 


scalar integrals 


analytical 


numerical 


Drell-Yan, DIS 


2 ^ 1 


3-loop, 3-point, 1 mass 


no 


yes 


to N^LO 


2^2 


2-loop, 4-point, 1 mass 


no 


yes (*) 




2^3 


1-loop, 5-point, 1 mass 


yes 




e+e" jjjjjjjjj 


1 ^ 3 


2-loop, 4-point, 1 mass 


no 


yes 


to NNLO 


1 4 


1-loop, 5-point, 1 mass 


yes 




PP, PP -> jj, jj, "f-f 


2^2 


2-loop, 4-point 


yes 


yes (*) 


to NNLO 


2^3 


1-loop, 5-point 


yes* 





Table 3: Knowledge on scalar integrals needed for various processes. In the column "analytical" , 
results existing in the literature to 0{e^) are listed; those marked with J are known to all orders in 
e. The column "numerical" shows where our method can improve the situation. The asterisk (*) 
indicates that more powerful numerical integrators than the ones used in the present work would be 
needed to deal with thresholds. 



In Table (3) we list some integrals entering the calculation of virtual corrections to phenomeno- 
logically relevant processes. We indicate for which ones analytical results exist in the literature and 
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where our numerical method improves the present situation. For the calculation of the NNLO QCD 
corrections of e^e~ — > jjj,jjj,j"tl two-loop box graphs with one massive external leg are needed. 
As has been demonstrated above these diagrams can be calculated numerically with our method. 
This is also true for 3~loop 3-point functions needed for a N^LO calculation of deep inelastic scatter- 
ing (DIS) and the Drell-Yan process^. Although our algorithm produces integrable functions which 
in principle always allow for a numerical evaluation, the presence of thresholds inside the integration 
region worsens the convergence properties of an integration routine. These cases are marked with 
an asterisk (*) in Table 1. It means that more efficient numerical integrators than the ones we used 
would be needed. 

After having calculated the purely virtual corrections of say a 2 ^ iV process involving L-loop 
integrals, the latter have to be combined with the 2 ^ + 1 corrections which include {L — 1)- 
loop integrals and where one has to integrate over the extra (unobserved) particle. This integration 
produces soft/collinear e-poles which require the knowlegde of the respective {L — l)~loop integrals 
to 0{e^), at least in the corresponding IR regions, since they have to be combined with these singular 
phase space integrals. However, most of the existing analytical results are known only to 0{e'^). In 
the context of our method, expansion of the results for the virtual integrals up to higher orders in e 
constitutes no principle problem, only the size of the FORTRAN routines will increase. 

In conclusion, we have presented a simple, constructive subtraction algorithm to deal with IR 
divergent multi-loop integrals. Although the iterated sector decomposition produces in general a 
relatively large number of subsector integrals, we demonstrated that by means of automatization 
our method allows for the numerical computation of highly nontrivial Feynman diagrams of two- 
and three-loop type. We pointed out that with our numerical approach, certain higher order QCD 
calculations can be tackled without waiting for further developments of analytical methods. 



®We note that this refers only to the calculation of the partonic cross section; for a complete analysis, the splitting 
functions also have to be known to this order. 
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